Finite density QCD with a canonical approach 
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We present a canonical method where the properties of QCD are directly obtained as a function of the baryon 
density p, rather than the chemical potential p. We apply this method to the determination of the phase 
diagram of four-flavor QCD. For a pion mass ~ 350 MeV, the first-order transition between the hadronic 
and the plasma phase gives rise to a co-existence region in the T-p plane, which we study in detail, including 
the associated interface tension. We obtain accurate results for systems containing up to 30 baryons and quark 
chemical potentials p up to 2T. Our T-p phase diagram agrees with the literature when ^ < 1. At larger chemical 
potential, we observe a "bending down" of the phase boundary. We compare the free energy in the confined and 
deconfined phase with predictions from a hadron resonance gas and from a free massless quark gas respectively. 



1. Introduction 

The last few years have seen remarkable 
progress in the numerical study of QCD at fi- 
nite chemical potential p. Still, the various meth- 
ods |1I2I3I4| suffer from systematic uncertainties, 
which limit their range of reliability to about 
^ < 1.0. For a recent review, see Ref. [5]. We 
try to address this apparent limitation by using a 
canonical approach |6I7| . where we focus on the 
matter density p, rather than the chemical poten- 
tial. Fixing the baryon number brings us closer 
to the experimental situation of a heavy-ion col- 
lision. Moreover, our method is particularly ap- 
propriate to explore few-nucleon systems at low 
temperature, and in principle, allows to study the 
bulk properties of nuclear matter and the nuclear 
interactions. Here, we extend its use and deter- 
mine the phase boundary between the confined 
phase and the quark gluon plasma, as illustrated 
in FigUJ 

The top figure shows a sketch of the conjec- 
tured phase diagram of QCD in the T-p plane. 
At small chemical potential, the phase transition 
is a rapid crossover, which ends in a second or- 
der endpoint, followed by a first order transition. 
Correspondingly, the bottom figure illustrates the 
transition in the T-p plane. The first order 
transition is manifested in a co-existence region. 
In this proceedings, we describe how we identify 
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Figure 1. Sketch of the conjectured QCD phase 
diagram in the grand-canonical (top) and canon- 
ical (bottom) formalism. 

the co-existence region and how we determine the 
phase diagram in the T-p as well as in the T-p 
plane. 
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2. Partition Functions 

We construct the canonical partition function 
Zc(T,Q) by fixing the number ofquarks ./V = 
J d 3 x ^(x) 70 ip(%) to Q. We insert a 5- 
function in the grand canonical partition function 
Z G c(T,p) 

J [DU\[D*][DV] e -s a lU;T]-s P [u,v,v;T] § ^ _ 

(1) 

The (5-function admits a Fourier representation 
6(N-Q} = fdp,j tfMX-Q) . We recognise 
ipi = ifliT as an imaginary chemical poten- 
tial and exploit the ^-^-periodicity [S] in pi of 
Zcc(T,p = ini) 

Z C (T,Q) = ^\\ dfne-^ZGci^mT) 



Q=3B 2ir 



-iZB- 



Z GC {T,m) (2) 



Thus, the canonical partition functions are the 
coefficients of the Fourier expansion in imaginary 
p of the grand canonical partition function. As 
a consequence of the ^^-periodicity, the canon- 
ical partition functions are zero for non-integer 
baryon number B = Q/3. 

From the canonical partition functions, the 
grand canonical partition function can be recon- 
structed using the fugacity expansion (in fact a 
Laplace transformation) 



Zac(T, p) 



dpe 3V »*Z c (T,p) 

OO 

dp e -Um P )-^P) (3) 



with the baryon density p = y and the Helmholtz 
free energy density f{T,p) = log Z C (T, p). 
The relation between baryon density and chemi- 
cal potential can be expressed as (p)(p) or p(p): 



1 



Zgc(T, p) 



or 



/ \ 1 df(j>) 
^ P) = 3^p- 



dppe 3V ^Z c (T,p) (4) 



(5) 



While the first expression is exact in any volume, 
the second is obtained via a saddle point approx- 
imation (exact in the thermodynamic limit) and 
may have more than one solution when solving for 
the baryon density at a given chemical potential, 
see FigO We discuss this issue in detail below. 

3. Method 

Following J]], we express the canonical parti- 
tion function in a ratio, which can be measured by 
Monte Carlo simulation as an expectation value: 

Zc(B,(3) 



Zgc{(3o = f3,H = i^i a ) 
7 { l . v f[DU] e- s ^Mdet(U^p Io ) 

X J_ f c -i3B^ det(£/;i/xj) 

27ri_, \T ) ' det{U;ip Io ) 

I— j(E±\ ~-i3BH*. det(U;ifn) . 
{ 2kJ_/\t) det(C/;w ) ;/3 °'^ 



, z c (U;B) 

"det(U;ip Ia )' Po ^ 10 



(6) 



where Zgc^Qi^I^Io) is the grand canonical par- 
tition function sampled by ordinary Monte Carlo 
methods, here for notational simplicity at /3q — (3. 
The Zc(U; B)'s are the Fourier coefficients of 
the fermion determinant for a given configura- 
tion {£/}. Although the average in Eq.@ should 
be real positive, the individual measurements are 
complex, with a sometimes negative real part. 
This is how the sign problem manifests itself in 
our approach. Moreover, a reliable estimate de- 
pends on a good overlap of our Monte Carlo en- 
semble with the canonical sector B at tempera- 
ture j3. We address this issue by following the idea 
of Ref. yP and including both confined and decon- 
fincd configurations in our ensemble. Indeed, we 
supplement the ensemble at (/3 c (a* = 0),/i = 0) 
with additional critical ensembles at imaginary 
chemical potential, non-zero isospin chemical po- 
tential, and ensembles generated with an asym- 
metric Dirac coupling 4 - in principle, any en- 
semble is allowed. We then combine all this 
information about a particular canonical parti- 
tion function Zg(B,[3) by Ferrenberg-Swendsen 
reweighting |1(J| . 
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The Fourier-coefficients of the determinant 
Zc(U;Q) are calculated exactly In the tem- 
poral gauge (C/4(x, t) = 1 except for £ = A^ — 1), 
the staggered fermion matrix M in the presence 
of a chemical potential can be written in the form 



Bo 
-1 




1 

-1 





1 

#2 







B 



P = 



N t -1 

n 

3=0 



where the B^s contain all space-like contribu- 
tions. Ref. m showed that the determinant can 
be computed for any chemical potential at the 
cost of diagonalising the so-called "reduced ma- 
trix" P. The determinant is given in terms of P's 
eigenvalues Ai , . . . , Aqv , where V is the spatial 
volume: 



6V 



det M(U; n) = e 3V ^ aNt J| (A; 



e 



Q=3V 

= Z c (U-Q)e-^ aN < . (7) 

Q=-3V 

Matching term by term, we then solve for the 
Fourier coefficients Zc(U;Q). This delicate step 
requires a special multi-precision library. The di- 
agonalisation of the reduced matrix P is compu- 
tationally intensive and takes 0(V 3 ) operations. 

4. Results 



this quantity with ^-rp- following Eq. JSJ. This 
assumption will be tested in FigO 
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Figure 2. (top) A sketch of the "scans" in our 
phase diagram, (bottom) The derivative of the 
free energy at fixed temperature as a function 
of the baryon number (baryon density). In the 
saddle point approximation, the y-axis is M; = 

F(B)-F(B-1) 
3T 



We study the Helmholtz free energy F(B) = 

Z ( B) 

— riog Zo(o) m a theory of four degenerate fla- 
vors of staggered quarks with mass ma — 0.05 



( l 



0.2, ~ 350 MeV) on a small, coarse 



6 3 x 4 lattice with volume ~ (1.8fm) 3 . For the 
quark mass we chose, the phase transition is first 
order at fi = 0, and presumably remains first or- 
der for all chemical potentials. 

We "scan" the phase diagram by varying the 
baryon density at fixed temperature, see Fig[21 
We measure F ^ B - > ' 3 ^ B ~ 1 - ) and assume the valid- 
ity of the saddle point approximation to equate 



Note that accurate results are obtained, up to 
high densities (> 5 baryons/ffn 3 ) and large chem- 
ical potentials (y ~ 2). The first order phase 
transition and the associated metastabilities are 
clearly visible in the "S-shape" of £ (B) . The low- 
density regime can be reasonably well described 
by a simple hadron resonance gas Ansatz, ^ = 
3/(T) sinh(3y ) with f(T) as the only free param- 
eter. The high-density regime almost corresponds 
to a gas of free massless quarks ^ = Nf 



N f 



v'i'/ 



f^) when taking cut-off corrections m t° 
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account. The solid line in Fig|2Ibottom) is ob- 
tained by fitting the linear and cubic terms in 
this expression. Instead of the free value 1, the 
fitted coefficients are 0.82(2) and 1.94(6) respec- 
tively. Thus, the equation of state for the quarks 
in the plasma phase differs little from the Stefan- 
Boltzmann law. This has been observed also in 
[H] and H2|. Such a result is surprising, given 
the experimental evidence for a strongly interact- 
ing plasma. Perhaps it simply reflects the lack of 
sensitivity of the quark equation of state to the 
coupling strength. Note in particular the striking 
independence of ^-(^) on temperature, as long 
as one is in the plasma phase. 

For a given temperature T, we identify the 
boundaries p\ and p 2 of the co-existence region 
and the critical chemical potential p c as follows. 
Equality of the free energy densities in the two 
phases, f(px) - 3p c pi = f(p 2 ) - 3p c p 2 , implies 

/ dp(f(p) - 3/ic) = . (8) 
J pi 

Since f'(p) is the quantity measured in Fig|21 
we determine p\ , p 2 and p c by a "Maxwell con- 
struction" illustrated in FigOU (top) for the tem- 
perature 3r — 0.92. The value of U; defining 
the horizontal line is adjusted to make the areas 
of the two "bumps" in the S-shape equal. The 
two outermost crossing points define p\ and p 2 , 
the boundaries of the co-existence region. Here, 
= 1.06(2) is the value of the critical chemical 
potential. 

We can cross-check this result by making use 
of the fugacity expansion Eq.0, see FigEJ (bot- 
tom) . For a given chemical potential, we measure 
the baryon number (B)(p). We see a rapid varia- 
tion at the same value ^ « 1.06, but the round- 
ing due to finite size effects is very strong. In 
contrast, our criterion for criticality (equality of 
the free energies) gives a clear first-order signal. 
In the thermodynamic limit, both approaches will 
indicate the same jump at p c . 

Thus, the S-shape we observe is a finite-size ef- 
fect. Its origin is clarified in Fig0] In a system 
of spatial size L, it takes a free energy 2L 2 a/T to 
create two planar, static interfaces separating the 
two phases, where a is the interface tension. Once 
these two non- interacting interfaces are present, 
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Figure 3. (top) The Maxwell construction allows 
to extract the critical chemical potential and the 
boundaries of the co-existence region, (bottom) 
Comparing the saddle point approximation (red) 
with the fugacity expansion (blue) . Strong finite- 
size effects in the latter obscure the first-order 
transition. 



the baryon density can be varied at constant free 
energy by changing their separation. Thus, the 
free energy required for the creation of these inter- 
faces is equal to the shaded area in Fig0J A crude 
trapezoidal integration gives a reasonable value 
for the reduced interface tension \fa]T ~ 35 — 45 
MeV. The contribution of these interfaces to the 
free energy density goes as 1/L, and vanishes in 
the thermodynamic limit. 

In Fig0we present the phase diagrams in the 
T-p as well as in the T-p plane. The top figure 
summarizes results from various methods, all for 
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rho/T 3 



line is to fall within the range a[i c = 0.35.. 0.43 
at (3 = 0, as predicted from a strong coupling 
analysis [T5] . 

FigEJbottom) shows the first phase diagram 
in the T-p plane. The densities at the bound- 
aries of the co-existence region seem to remain 
constant for T < 0.85T C already, with pqgp = 
1.8(3)S/fm 3 and p C onfincd = 0.50(5)5/fm 3 . The 
latter is a plausible value for the maximum nu- 
clear density in our 4-flavor, — 350 MeV QCD 
theory. 



Figure 4. Connection between the S-shape of the 
derivative of the free energy and the interface ten- 
sion: the shaded area represents the free energy 
required to create two planar interfaces. 



the same theory: 4 flavors of staggered quarks 
with am = 0.05, N t = 4 time-slices; only the spa- 
tial volume varies as indicated. All methods agree 
nicely up to fi/T ~ 1, but have different systcm- 
atics which become dominant at larger chemical 
potentials. The multi-parameter approach of PP 
uses a single Monte Carlo ensemble generated at 
(PcP = 0), and determines the phase transition 
from Lee- Yang zeroes. To understand the system- 
atic errors involved, we have repeated this study, 
but used instead the peak of the specific heat to 
determine the phase transition. Our results are 
completely consistent with U|. However, we have 
noticed that the "sign problem" grows dramati- 
cally more severe with increasing chemical poten- 
tial, but that our statistical error, based on jack- 
knife bins as in [Q, fails to reveal this breakdown 
of the method. 

The parabolic fit [Hj is consistent with the more 
elaborate approach of 0] . Both methods perform 
an analytic continuation from imaginary /i, for 
which the systematic errors are hard to quantify. 
Our new results are shown in red. There is no 
strong inconsistency with other results, but we 
observe a clear sign of bending down starting at 
^ ~ 1.3. In fact this must happen, if the critical 
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Figure 5. (top) The phase diagram in the T-\x- 
plane, obtained by 3 different methods. The two 
arrows indicate the range of the strong coupling 
((3 — 0) prediction \T^\. (bottom) The phase dia- 
gram which we obtain, in the T-p-plane. 
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5. Conclusions 

We study QCD in a canonical framework, 
which is promising for the study of few-nucleon 
systems at low temperature, but proves also capa- 
ble of exploring high density regimes (/i/T < 2) at 
temperatures T > 0.8T C . We have determined the 
phase boundary between the confined phase and 
the quark gluon plasma in both the T-p and the 
T-fi plane. In the latter, our results are in agree- 
ment with the literature, however we observe a 
bending down of the critical line at ^ ~ 1.3. The 
two phases can be rather well described by the 
hadron resonance gas at low densities and by a 
weakly interacting massless gas at high densities. 
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